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Unraveling the interplay between connectivity and spatio-temporal dynamics in neuronal net¬ 
works is a key step to advance our understanding of neuronal information processing. Here we 
investigate how particular features of network connectivity underpin the propensity of neural net¬ 
works to generate slow-switching assembly (SSA) dynamics, i.e., sustained epochs of increased firing 
within assemblies of neurons which transition slowly between different assemblies throughout the 
network. We show that the emergence of SSA activity is linked to spectral properties of the asym¬ 
metric synaptic weight matrix. In particular, the leading eigenvalues that dictate the slow dynamics 
exhibit a gap with respect to the bulk of the spectrum, and the associated Schur vectors exhibit 
a measure of block-localization on groups of neurons, thus resulting in coherent dynamical activ¬ 
ity on those groups. Through simple rate models, we gain analytical understanding of the origin 
and importance of the spectral gap, and use these insights to develop new network topologies with 
alternative connectivity paradigms which also display SSA activity. Specifically, SSA dynamics in¬ 
volving excitatory and inhibitory neurons can be achieved by modifying the connectivity patterns 
between both types of neurons. We also show that SSA activity can occur at multiple timescales 
reflecting a hierarchy in the connectivity, and demonstrate the emergence of SSA in small-world 
like networks. Our work provides a step towards understanding how network structure (uncovered 
through advancements in neuroanatomy and connectomics) can impact on spatio-temporal neural 
activity and constrain the resulting dynamics. 
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Neural networks display a wide range of spatio- 
temporal behaviors. Understanding how this complex 
orchestration of neuronal firing activity is determined by 
the structure of the network ( i.e ., its wiring) is an im¬ 
portant step towards comprehending how neural com¬ 
putation is manifested, especially given the growing ex¬ 
perimental access to temporal recordings and connec¬ 
tomics. Here we investigate the link between network 
structure and the dynamics of neuronal assemblies in the 
context of leaky-integrate-and-fire (LIF) networks. We 
show how structural features in the wiring of the net¬ 
work can introduce additional time-scales to the dynam¬ 
ics, and how such structured wiring can lead to spatio- 
temporally segregated, coherent activity of groups of neu¬ 
rons. Using rate models we gain insight into how such 
spatio-temporal dynamics emerge as a direct consequence 
of the spectral properties of the network, and use this 
understanding to examine further circuit topologies that 
enable phenomenologically similar behavior, yet rely on 
fundamentally different connectivity and functional pat¬ 
terns. 
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I. INTRODUCTION 

Neuronal ensembles exhibit a broad repertoire of ac¬ 
tivity patterns. Such dynamics are governed by a time- 
evolving network of synaptic connections with an intri¬ 
cate, yet structured, organization. Due to the advance¬ 
ment of connectomics, our knowledge about such net¬ 
works is rapidly growing, and increasingly detailed maps 
of neuronal wiring are becoming available. In parallel, 
modern recording techniques, such as calcium imaging 
and multi-electrode arrays, allow neuroscientists to mon¬ 
itor the activity from thousands of neurons simultane¬ 
ously, with recordings from entire brains at single neuron 
resolution becoming technologically feasible [1-3]. The 
observed dynamics of neural networks exhibit an inter¬ 
play of structured spatio-temporal scales, which underpin 
a wide range of cognitive functions [4]. 

The idea that neuronal group activity induced by net¬ 
work structure is at the core of neural computation dates 
back at least to the work of Hebb [5], who hypothesized 
that the transient activity of groups of neurons (so called 
cell assemblies ) is the currency of information processing 
[4, 6]. This notion is supported by recent experiments 
showing that reciprocal connections between neurons oc¬ 
cur above chance level [7, 8], especially if neurons receive 
common inputs [9, 10]. In the case of the visual system, 
for instance, excitatory neurons with similar response fea¬ 
tures tend to be more connected to each other [11, 12]. 
Moreover, studies have demonstrated that neurons ex¬ 
hibit layer-specific connectivities within rodent sensory 
cortex [13] and neocortex [14]. In addition, organized ar¬ 
chitectures have been observed to occur at multiple hier- 
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archical scales [15-18] and in non-mammalian organisms 
[19]. These findings suggest that cortical regions contain 
well-defined subnetworks. However, the underlying ques¬ 
tion is whether given a network topology, we can predict 
the potential of the network to sustain structured spatio- 
temporal activity. Such questions are not only of inter¬ 
est for network dynamics, but have also implications for 
memory formation and learning, since neural networks 
undergo topological changes over time due to plasticity 
[ 20 , 21 ]. 

Recently, it has been shown computationally (see 
e.g. Ref. [22]) that leaky-integrate-and-fire (LIF) net¬ 
works with equal excitatory and inhibitory connection 
net strengths (i.e. balanced [23]) yet with clustered ex¬ 
citatory connections, can exhibit prolonged heightened 
group activity, with the activity transitioning between 
groups in the network (Fig. 1A). Here we characterize the 
emergence of such slow-switching segregated dynamics in 
balanced LIF networks as a result of the network connec¬ 
tivity. Specifically, we find that the spectral properties 
of the synaptic weight matrix (i.e., the existence of an 
eigenvalue gap and a block-localized dominant subspace) 
provide a criterion to predict the appearance of such ac¬ 
tivity in the network. We then use simple linear rate 
models to gain insight into the mechanisms underpin¬ 
ning the origin of such dynamics in structurally clustered 
LIF networks. Using these insights, we construct novel 
LIF topologies that display slow-switching group activity 
with distinct properties: involving both inhibitory and 
excitatory neurons; exhibiting multiple slow time-scales; 
as well as demonstrating the possibility of such dynam¬ 
ics in networks with no obvious clustered connectivity, 
such as small-worlds. Finally, we discuss briefly possible 
implications of the different wiring schemes for neural 
computation. 


II. RESULTS 

A. Slow-switching assemblies in LIF networks with 
clustered excitatory connections: spectral insights 

Clustered excitatory topologies in a balanced LIF net¬ 
work can lead to dynamics in which localized high activ¬ 
ity states transition between assemblies of neurons within 
the network [22]. This is illustrated in Figure 1A, where 
the dynamics of an unstructured and a balanced clus¬ 
tered network with 20 groups are shown side by side (see 
Materials and Methods for a description of the networks). 
Hereafter, we will refer to such activity as slow-switching 
assembly (SSA) dynamics. Visually, SSA dynamics man¬ 
ifests itself as bands of increased activity in the raster 
plots, and can be statistically quantified from the result¬ 
ing spike-train dynamics a posteriori (see Eq. (18) in Ma¬ 
terials and Methods). Ideally, however, we would like to 
establish a priori , solely from the given connectivity, the 
possibility of such dynamical patterns emerging. 

The full dynamics of LIF networks are notoriously 
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Figure 1. Network dynamics and eigenvalue spectra of two 
LIF networks: one with uniform synaptic connections (left), 
and one with 20 groups of clustered excitatory connections 
(right). To create the clustered network, excitatory neurons 
were partitioned into groups (with in-group connection proba¬ 
bility p^ and out-group connection probability p < p ^ E ) 
while keeping the average connectivity constant (see Materials 
and Methods and Ref. [22]). The ratio Ree — PtL E / Pout con¬ 
trols the strength of the excitatory clustering. A Visualization 
of the network topologies (top) and exemplars of raster plots 
(bottom). The dynamics of the clustered network exhibit the 
banded structure associated with slow-switching group activ¬ 
ity. The magnitude of the activity can be characterized statis¬ 
tically a posteriori from the data through the spike-rate vari¬ 
ability metric S, defined in Materials and Methods Eq. (18), 
as discussed in the text. In this case, the unclustered network 
has S = 0.035 while the clustered network has a much larger 
value S = 8.23. B Eigenvalue spectra of the network weight 
matrices W. The weighted connectivity matrix of the clus¬ 
tered network exhibits a clear eigengap AA separating the 19 
eigenvalues with largest real parts from the cloud of eigenval¬ 
ues in the bulk. There is no such eigengap for the unclustered 
network. As indicated by the two arrows, both matrices have 
a pair of complex conjugate eigenvalues associated with the 
(damped) global activation modes of the networks character¬ 
istic of balanced networks (see text and Ref. [24]). 


difficult to analyze due to their inherent non-linearity; 
hence an exact analytical treatment of the dynamical 
evolution for an arbitrary clustered topology is essen¬ 
tially intractable. However, two concepts from spec¬ 
tral graph theory and linear systems provide valuable in¬ 
sights: (i) for symmetric, non-negative connectivity ma- 
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trices, it can be shown that a modular network structure 
implies a gap in the spectrum of the graph (i.e., in the 
set of eigenvalues of the weight matrix) [25, 26], as well 
as the block-localization of the associated eigenvectors on 
the modules of the network (noting that isolated eigenval¬ 
ues may also be the result of other features [27]); (ii) for 
linear systems, a gap in the spectrum of the graph results 
in a separation of time scales in the dynamical process 
[24, 28, 29]. This relation between the modular struc¬ 
ture, the eigenvalues and associated eigenvectors, and 
linear network dynamics can be used to discover mod¬ 
ular structures in graphs from a dynamical perspective 
[30-33]. 

In fact, the weighted connectivity matrices of unclus¬ 
tered and clustered LIF networks [22] display different 
spectral characteristics. Figure IB shows the spectra of 
two networks with 1600 excitatory and 400 inhibitory 
neurons with unclustered (left) and clustered (right) 
topologies. In the unclustered case, we find the expected 
circular distribution of eigenvalues, which follows from 
the properties of random graphs [34, 35], although the 
presence of groups of neurons with different cardinalities 
and variances means that the eigenvalue distribution is 
not completely uniform on the circle [34]. Note also that 
the balanced construction of the LIF network, with a 
marginally larger inhibitory input for each neuron in or¬ 
der to keep the network stable, leads to the existence of 
one pair of complex conjugate eigenvalues which lies sep¬ 
arate from the main bulk of the spectrum (black arrows 
in Figure IB). This eigenvalue pair is associated with 
the global activation mode of the network (as explained 
below in Section IIB and in Ref. [24]). As shown in Fig¬ 
ure 1A, this unclustered network exhibits the expected 
asynchronous, unstructured neuronal spiking dynamics. 

In contrast, the LIF network with clustered excitatory 
neurons displays banded SSA dynamics. Spectrally, its 
weight matrix exhibits a clear gap A A along the real axis 
of its spectrum (Figure IB, right) and, as shown below in 
Section IIB 3 (Fig. 5), the associated Schur vectors also 
exhibit a measure of structural block-localization. We 
remark that, in general, this should not be expected a 
priori , since the weight matrices are asymmetric and in¬ 
clude both positive (excitatory) and negative (inhibitory) 
couplings. 

To ascertain the dynamical relevance of the spectral 
gap, we examined the relation between the clustering 
strength in the network, defined as the ratio of proba¬ 
bilities of connections inside and outside the neural as¬ 
semblies ( Ree = Phi E /PoJt )> the spectral gap (AA), and 
the magnitude of the numerically observed SSA dynam¬ 
ics. To quantify the assembly spike-rate variability, we 
have defined two complementary metrics. First, the met¬ 
ric S measures the heterogeneity in the firing rates of the 
putative cell-assemblies averaged over the simulation (see 
Eq. (18) in Materials and Methods). A large value of S 
indicates that the average firing rates of the assemblies 
are diverse, whereas a low S indicates that all groups 
have very similar firing rates at all times and no group 




Ree 

0.7 

0.6 

0.5 

AA o.4 

0.3 

0.2 

0.1 


10 

8 

6 

4 

2 


l 


0 0.1 0.2 0.3 0.4 0.5 

AA 


i 11 1 » i i 11 « » i ■ 


2 2.5 3 3.5 4 4.5 5 


Ree 


Figure 2. The relationship of observed SSA dynamics with 
the structural connectivity clustering of the LIF network, 
Ree , and the spectral gap, AA. The presence of SSA dy¬ 
namics is quantified through the spike-rate variability metric 
(S'), which measures the variance of the firing rates of the 
assemblies normalized by a randomly shuffled bootstrap: S 
increases with increasing SSA activity and S ~ 0 for com¬ 
pletely asynchronous activity (as in the unclustered case in 
Figure 1A) (see Materials and Methods). A Spike-rate vari¬ 
ability is plotted as a function of Ree (left) and A A (right) for 
different network sizes (dots: raw data from simulations, line: 
mean, shading: standard deviation). Above a certain cluster¬ 
ing threshold, SSA emerges and increases as Ree grows; the 
intensity of the SSA dynamics is in line with the presence of 
an eigenvalue gap AA in the weight matrix. B Spike-rate vari¬ 
ability as a function of Ree (left) and A A (right) for a network 
of 2000 neurons with different numbers of clusters, yielding 
qualitatively similar results (dots: raw data from simulations, 
line: mean, shading: standard deviation). C Relationship be¬ 
tween the clustering strength Ree and the spectral gap A A. 
Observe that Ree is not sufficient to determine AA, i.e. AA 
is influenced by other aspects such as the network size and 
number of groups. 


shows elevated firing. Under SSA dynamics, the variabil¬ 
ity of firing rates across groups increases in time as the 
assemblies transition between high and low firing rates. 
As discussed below in Section IIB 4, it is possible that 
the heightened firing activity is localized in a particular 
assembly and does not switch between assemblies. This 
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Figure 3. Structure of the eigenspectra of stylized and real¬ 
istic networks that can exhibit SSA activity. A Schematic of 
the stylized model of a clustered network (l)-(2) and illus¬ 
tration of the corresponding eigenvalues (6). B Schematic of 
a more realistic clustered LIF network with the same wiring 
scheme on many groups and multiple neurons per group with 
an illustration of the associated spectrum. 


non-switching dynamics where only a group of neurons 
exhibits elevated firing over the entire simulation time 
can be detected using the second variability measure St , 
defined in Eq. (20). For SSA dynamics, both S and St 
give similar results (see Fig. SI). In the rest of the pa¬ 
per, we focus on SSA dynamics and mainly use S', except 
when discussing the transition to non-switching behavior 
in Section IIB 4. Definitions of these quantities are given 
in Materials and Methods. 

Figure 2 shows that SSA becomes observable above a 
threshold of the clustering strength Ree , although the 
dependence of this threshold on the size of the network 
and number of clusters does not seem to follow an obvious 
pattern. On the other hand, our numerics show that the 
spectral gap provides a more direct indicator of the pres¬ 
ence of SSA dynamics in the network. The relationship 
between Ree and AA (Figure 2C) shows that knowing 
Ree is not sufficient to determine A A as the spectral gap 
depends on other factors such as the network size and the 
number of groups in the network. Together with the ex¬ 
amination of the structure of the associated orthogonal 
subspace (see Section IIB 3), such spectral characteriza¬ 
tion may be used to establish the potential of networks 
to sustain SSA activity, even if there is no obviously clus¬ 
tered network model known a priori. 

The observed spectral properties of clustered LIF net¬ 
works lead us to consider a stylized linear rate model for 
neuronal activity in the next section, as a means to gain 
insights into the defining factors of the wiring diagram 
leading to SSA. 


B. Linear rate models and slow localized activity 
in clustered LIF networks 

1. A stylized linear rate model for networks with clustered 
excitatory neurons 

Motivated by the above findings, we proceeded to in¬ 
vestigate how much of the observed LIF dynamics can 
be described by simple rate models. From above, we ex¬ 
pect that the spectral properties of the network lead to 
SSA dynamics that is coherent over, and transitioning 
between, clusters. Hence we consider simple rate models 
on coarse-grained networks where each node describes 
the behavior of a cluster. As will become clear below, 
this reduced description allows us to identify the mech¬ 
anism underpinning the onset of SSA activity without 
having to consider a plethora of statistical and technical 
detail. Furthermore, our main conclusions can be eas¬ 
ily extended to multiple groups of clustered excitatory 
neurons, or to scenarios with probabilistic connectivity 
between the nodes (see e.g. Ref. [24]). 

To explore this idea in the simplest setting, consider 
a linear firing rate model with three groups: two groups 
of clustered excitatory neurons and a group of inhibitory 
neurons with uniform coupling. This wiring scenario is 
abstracted in the form of a three-node network (Figure 
3A), where each node represents a group of neurons, and 
a simple rate model is given by the following equations: 

dv 

r— = -r + Wr + £ = -(/ - W )r + £. (1) 

Here r = (r e i, r e 2 , r^) T are the firing rates of the neuron 
groups; W is the synaptic connectivity weight matrix; 
and £ is a random input to the network. The firing rates 
are measured relative to some baseline activity, and may 
be positive or negative. 

To describe this topology with two excitatory clusters 
we define a connectivity matrix W of the form: 


1 s e 

—kw\ 


W = [ e s 

—kw I , 

( 2 ) 

\w/ 2 w/2 

—kw J 


s > e > 0, k > 1, 

w = (s + e). 

( 3 ) 


The first inequality in (3) guarantees that the coupling 
strength is positive: s — e > 0, i.e., the cross-coupling 
between the excitatory groups is weaker than the intra¬ 
group coupling, as would be expected from the notion of 
a cluster. The second inequality in (3) ensures a balanced 
network (while keeping k ~ 1), such that each group has 
at least the same amount of inhibitory and excitatory 
input: Wij = (1 — k)w < 0, Vi. 

To understand the dynamical behavior of this system, 
we evaluate the spectral properties of the connectivity 
matrix W [24, 29, 36]. Note that we are dealing with 
a non-normal system, in which the eigenvectors may not 
be orthogonal and thus will provide a possibly misleading 
description of the dynamics [37]. We therefore consider 
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the Schur decomposition of W, which yields an orthonor¬ 
mal basis of the system [24, 36], i.e., it finds an orthog¬ 
onal matrix U and an upper triangular matrix Q such 
that W = UQU T . Here U contains an orthonormal ba¬ 
sis while Q contains the eigenvalues of W on the diagonal 
and non-zero elements above the diagonal accounting for 
effective feedforward connections between the different 
modes [24, 36]. 

For the matrix (2), the associated upper triangular 
Schur form Q is: 

( —Wff 0 \ 

o o J, (4) 

w + = w(k — 1 ), Wff = (fc + 1/2 )(w + e), 

and the matrix containing the orthonormal (Schur) basis 
is U = (ui U -2 U 3 ) given by 



( 5 ) 


Each of the Schur vectors Ui represents a pattern of rates, 
a particular mode of firing on the network. The first 
mode ui is associated with an overall mean firing pat¬ 
tern, while the second mode 112 corresponds to a rela¬ 
tive difference in firing between inhibitory and excita¬ 
tory groups. Note that in both cases the firing patterns 
of the two groups of excitatory neurons are aligned. In 
contrast, the third mode 113 describes an antagonistic 
activity localized on the excitatory neuron groups (i.e., 
when the firing rate of one excitatory group increases the 
other decreases, and vice versa) with the inhibitory node 
remaining unaffected at a baseline firing rate, precisely 
in line with the SSA behavior observed in the full LIF 
network (Fig. 1A). Observe also that, while the first two 
modes are coupled via an off-diagonal term in Q (i.e., an 
effective feedforward connection [ 24 , 36 ]), the switching 
behavior described by the third mode is uncoupled from 
the other system modes. 

The structure of the Schur decomposition provides 
us with insight into the dynamics of the model (1). 
From ( 4 ), the eigenvalues of W are: 

Ai = — w + < A 2 = 0 < A 3 = s — e, ( 6 ) 

whence the eigenvalues governing (1) are 771 = —1 — re + , 
772 = — 1 and 773 = — 1 + (s — e). For the system to be 
linearly stable, we require 773 < 0, which implies that the 
clustering strength in a stable system is constrained to 
be 

0 < s — e < 1. (7) 

Without an external input, the three modes decay expo¬ 
nentially with time constants r* = —l/rji = 1/(1 — A i). 
Therefore, it is possible to slow down the time scale asso¬ 
ciated with the SSA mode 113 by increasing the cluster¬ 
ing strength s — e 1 . The perturbations introduced by 
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Figure 4. SSA dynamics can be achieved by changing solely 
the synaptic strengths even when the topological connections 
are kept uniform. To create clustered networks, the ratio 
Wee = / Wont was varied, where w^ E and w E f^ refer 

to the in-group and out-group synaptic weights, respectively. 
The spike-rate variability S measuring the intensity of SSA 
dynamics is shown as a function of Wee (left - dots: raw data 
from simulations, line: mean, shading: standard deviation) 
and the spectral gap A A (right). Note that the connection 
probability is uniform for all the simulations, i.e. Ree = 1 , 
and only the clustering of the weights is varied. 


the random input induce mode-mixing and, in particular, 
break the symmetry between the firing rates of the two 
groups of excitatory neurons. These perturbations ex¬ 
cite the SSA mode u 3 , which will decay with time scale 
1/(1 — A 3 ) towards the solution with uniform rates on all 
groups (r* = 0 ). 

In conclusion, the stronger the clustering (i.e., the 
larger s — e < 1 ), the slower the decrease of the asym¬ 
metric transients, and the more prominent SSA dynam¬ 
ics becomes. Note that this slow time scale in the ac¬ 
tivity patterns is directly associated with the eigengap 
A 3 — A 2 = s — e (Fig. 3), which is dictated by the clus¬ 
tering strength, in agreement with our observations re¬ 
lating the presence of SSA with AA in Figure 2 . The 
insights gained from this simple model have testable im¬ 
plications for LIF networks and lead to further ideas for 
neuro-physiological network wiring structures, which we 
explore in the following sections. 


2. The linear rate model and the full dynamics of clustered 
LIF networks 

Our rate model can be extended to networks with c 
groups and the overall structure of the eigenspectrum of 
W will essentially retain its features (Figure 3; see [24] 
for a related discussion). First, there will be a ‘nega¬ 
tive’ eigenvalue (Ai in ( 6 ) related to the pair Aeaiance in 
Fig. 3B), which reflects the fact that the network is bal¬ 
anced and stable, i.e. the associated global activity mode 
decays. Second, there will be a bulk of ‘small’ eigenvalues 
centered around the origin (A 2 in ( 6 ) and in general the 
set {Aeuik} in Fig. 3B), which stem from the random con¬ 
nectivity present in the network—a random matrix has 
a circular eigenvalue distribution around the origin). Fi¬ 
nally, there will be an eigenvalue gap separating a small 
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set of c — 1 ‘large positive’ eigenvalues (A 3 in ( 6 ) and 
{^cluster} in Fig- 3B) from the bulk of the spectrum. 
These are the ‘slow’ eigenvalues associated with an or¬ 
thogonal Schur subspace with a block structure coherent 
with the clusters, and are thus responsible for the ap¬ 
pearance of a new time scale in the system originating 
the observed SSA dynamics. This structure of the eigen¬ 
values can be seen in the spectrum of the LIF network 
in Figure IB. 

Inspired by the analysis of the linear rate model, we 
investigate how the results translate to the full nonlin¬ 
ear LIF dynamics. First, it is clear from above that only 
the difference between the effective intra- and inter-group 
coupling strengths is essential for the observed slow activ¬ 
ity. Such difference can be the result of clustered topo¬ 
logical connectivity as in Ref. [ 22 ], but can be equally 
achieved keeping the topology uniformly connected and 
tuning the synaptic weights. To assess this possibility, 
we conducted numerical simulations of LIF networks with 
uniform connection probabilities between excitatory neu¬ 
rons, yet with larger intra-group synaptic weights while 
keeping the average weight constant (see Materials and 
Methods). Figure 4 shows that tuning the weights to¬ 
wards a clustered structure (as given by the weight ratio 
Wee — w hi E / w oJ!t ) also leads to the emergence of SSA 
dynamics linked to a spectral gap. Although unsurpris¬ 
ing from a linear systems perspective, this result confirms 
the applicability of the spectral characterization to non¬ 
linear LIF networks and, specifically, establishes its rele¬ 
vance for SSA dynamics in networks where only synaptic 
weights can be tuned. 


3. The block-localization of the dominant linear subspace 

We now consider a critical factor in the emergence of 
SSA, namely, that the leading Schur vectors of the weight 
matrix of LIF networks exhibit a measure of consistent 
block-localization on the groups of neurons which then 
exhibit cell assembly dynamics (Figure 5). Due to the 
non-normality outlined above [24, 36], we consider or¬ 
thogonal Schur vectors and not eigenvectors. A typical 
Schur vector corresponding to one of the eigenvalues in 
the leading group {Aciuster} exhibits a coherent pattern 
on each of the groups of the network, such that all neu¬ 
rons within a particular group are driven uniformly. This 
block-uniformity results in the grouped dynamics that 
characterizes SSA. Representative patterns are shown in 
Figure 5 A. In contrast, the Schur vectors from the bulk 
of the spectrum {Aeuik} show no specific pattern of lo¬ 
calization on any group of neurons. Note that the Schur 
vectors show no localization pattern on the inhibitory 
neurons either, in line with our analysis. The block- 
localization of the leading Schur vectors constitutes an 
additional spectral characterization to assess the emer¬ 
gence of SSA dynamics from the weight matrix alone. 

In Figure 5, we quantify how the patterns observed in 
the SSA dynamics align with the dominant Schur vectors 


of the weight matrix. To do this, we first performed a 
principal component analysis (PC A) of the simulated fir¬ 
ing rates of the full LIF network and extracted the top ac¬ 
tivation patterns {pi,..., p c } corresponding to the first 
principal components (see Figure 5B and Materials and 
Methods). To measure the alignment of the data patterns 
{pi} with the dominant Schur vectors {ui,..., u c } of the 
weight matrix W corresponding to the leading eigenval¬ 
ues above the gap, we computed the first principal angle 
6 between these two subspaces: the smaller the angle, the 
closer the alignment [38, 39]. Figure 5D shows that the 
dominant subspace of the Schur vectors of the weight ma¬ 
trix (W) are indeed aligned with the observed patterns 
and consistent with the embedded groups of neurons in 
the network. Finally, Figure 5E shows again that tuning 
the synaptic weights to cluster the network while keep¬ 
ing the topology uniformly connected renders a similar 
outcome. 


4. Increasing the clustering beyond the linearly stable 
regime: one dominant assembly 

Until now, our linear analysis has concentrated on the 
relevant regime for SSA, where the system is linearly sta¬ 
ble (7): A 3 = s — e < 1. However, it is possible to increase 
the clustering strength (s — e), so that the intra-cluster 
connections are much stronger than the inter-cluster con¬ 
nections, and the linear system has an unstable eigen¬ 
value. For the LIF network, if we increase Ree 01 * Wee 
so that A max G {Aciuster} > 1, then the associated local¬ 
ized firing mode, once excited, will activate itself. This 
regime can thus lead to a ‘winner takes all’ solution, 
where one cell assembly dominates the firing of the whole 
network and suppresses all other neurons [40]. 

From the perspective of the observed dynamics of the 
LIF network, an increase of the clustering strength ( Ree 
or Wee) initially leads to the development of the eigen- 
gap for {Aciuster} and, the emergence of SSA dynamics. 
However, when the clustering strength becomes too large 
(and the largest eigenvalue goes beyond 1 ), a single as¬ 
sembly starts to dominate the firing pattern and the fir¬ 
ing variability is reduced across time, an effect which is 
captured by a strong reduction in St, as illustrated in 
Figure 6 . This behavior, which has been observed pre¬ 
viously [ 22 ], is thus also closely linked to the spectral 
properties of the system. Note that the linear condition 
Amax > 1 is only indicative: the non-linearity and bound¬ 
edness of LIF dynamics can lead to (non-linearly stable) 
SSA dynamics beyond this condition. 


C. Beyond clustered excitatory neurons: SSA 
dynamics involving both excitatory and inhibitory 
neurons. 

So far we showed that a linear rate model can provide 
valuable insights into the mechanisms underpinning SSA 
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Figure 5. Unraveling the association between the leading Schur vectors of the weight matrix and the observed SSA activity. 
A Illustration of a weight matrix of the N — 2000 neuron LIF network with c = 20 groups of excitatory neurons (Ree = 3.4) 
shown together with the heatmap of the real parts of its first 25 Schur vectors. The first 19 Schur vectors exhibit block-uniform 
patterns relatively constant within each of the 20 groups of neurons, whereas no such pattern is observed for the other Schur 
vectors. Note that there is no pattern discernible over the inhibitory neurons. The leading Schur vectors correspond to the 
cloud of ‘slow’ eigenvalues above the gap, as indicated in C, and span a patterned dominant subspace that induces grouped 
dynamics in the network. B A long simulation (80s) of the LIF network dynamics (only the first 2s shown) was analyzed using 
PCA and the first 25 principal components (PCs) are shown. Reflecting the banded structure of the simulated dynamics, the 
leading PCs also show a block-patterned structure consistent with the neuronal groups. C On the spectrum of W , we indicate 
the group of leading eigenvalues above the gap associated with the dominant subspace. D The alignment between the dominant 
Schur subspace of the W matrix and the subspace of the strongest principal components is measured by the first principal angle 
0 (21). Above a threshold of the clustering strength Ree , both subspaces become highly aligned in line with the observations 
in Fig. 2 (dots: raw data from simulations; line: mean; shading: standard deviation). E The same effect is observed when the 
clustering is introduced in the weights by varying Wee , as in Fig. 4. 


dynamics in networks with clustered excitatory neurons. 
The crucial point for the emergence of SSA is the separa¬ 
tion of time scales, dictated by the splitting of the leading 
eigenvalues of the weight matrix, together with the block- 
localization of the associated dominant subspace. These 
spectral properties can be introduced into LIF networks 
by entirely different synaptic couplings, and we now con¬ 
sider a mechanism in which the inhibitory neurons are 
involved more centrally in generating SSA dynamics. 

Consider the wiring diagrams presented in Figure 7A- 
B, representing the simple rate model and its correspond¬ 
ing LIF schematic. In this case, the wiring does not 
involve clustered coupling between groups of excitatory 
neurons, but rather relies on preferential connectivity 
patterns between inhibitory and excitatory neuron groups 


only. Each group of excitatory neurons activates prefer¬ 
entially an associated group of inhibitory neurons and, 
in turn, this group of inhibitory neurons feeds back more 
weakly to its associated group of excitatory neurons (see 
Materials and Methods for a full description of the net¬ 
work). Therefore the effective functional circuitry con¬ 
sists of both excitatory and inhibitory neurons embedded 
in a feedback loop and, as a consequence, the inhibitory 
neurons play an integral role in generating the spatio- 
temporal dynamics and display SSA dynamics too, as we 
show below. 

This wiring mechanism was suggested by the styl¬ 
ized firing rate model with two coupled pairs of in¬ 
hibitory/excitatory feedback loops in Figure 7A. This 
system is described by Eq. (1) with a vector of firing 
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Figure 6. Effect of large clustering strength on SSA activity. 
A Large values of the clustering Ree lead to linear instability 
of the SSA dynamics and localization of the activity on one 
assembly. As measured by the spike rate variability across 
time (St), the increase of Ree leads to SSA (signalled by the 
increased value of St)- If the clustering increases further, St 
decreases, as the dynamics becomes dominated by one assem¬ 
bly only. The dots denote raw data from simulations (line: 
mean; shading: standard deviation). Inset: examples of raster 
plots for three data points in the three regimes. The analy¬ 
sis corresponds to a clustered LIF network of 1000 neurons. 
B Plot of the eigenvalue with the largest real component 
Amax against the clustering strength Ree - The linear condi¬ 
tion A max > 1 is a good indicator of the dynamics becoming 
dominated by one cell assembly. 


rates for the four groups r = (r e i, r e 2 , r* 2 ) T and a 
synaptic weight matrix defined as: 
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Similarly to (2), the system is balanced and s — e captures 
the clustering strength within the excitatory-inhibitory 
pairs. The Schur decomposition of W leads to the fol- 
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Figure 7. SSA dynamics can result from excitatory-to- 
inhibitory feedback loops. A Schematic of a stylized network 
model with two excitatory-to-inhibitory feedback loops corre¬ 
sponding to the model (8). B Illustration of the correspond¬ 
ing LIF network with such a co-clustered feedback mechanism 
between neuron types. C Raster plot of the dynamics of a 
co-clustered LIF network with N — 2000 neurons and c = 20 
pairs with excitatory-to-inhibitory feedback. Note that the 
inhibitory neurons also exhibit SSA dynamics here. D Spec¬ 
trum of the weight matrix exhibiting a spectral gap both on 
the left and right hand side of the bulk of the spectrum (see 
text for details). 
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where the eigenvalues of W are on the diagonal. When 
the leak term is considered, the largest eigenvalue of ( 1 ) 
is (—1 + Vk(s — e)). Hence, to keep the system stable we 
need to constrain 


0 < Vk (s — e) < 1 , 


( 11 ) 


and the spectral gap is again controlled by s — e. In the 
associated orthonormal Schur basis, the first two modes 
of the firing rate dynamics 


ui = 



u 2 = 



( 12 ) 


are again global ‘sum’ and ‘difference’ modes, which in¬ 
teract via a balanced amplification mechanism [24]. How¬ 
ever, there are also two localized modes 113 and 114 that 
can lead to slow structured activity: 


1 


u 3 = 


yj2k + 2 


fVk\ 

— yfk 
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V -1 / 


u 4 = 


\/ 2 k 2 


1 \ 

-1 

-Vk 

WkJ 


(13) 

Of these, u 3 is associated with the largest eigenvalue 
and describes the slow(est) dynamics. This mode cor¬ 
responds to a firing pattern of correlated activity within 
the pairs (r e i,rn) and (^ 2 ,^ 2 ), and anti-correlated ac¬ 
tivity across the pairs. As before, this analysis extends 
to networks with more groups and/or stochastic coupling 
between groups (for related arguments see supplementary 
material of Ref. [24]). 

In Figure 7C-D, we show that the implementation of 
this wiring mechanism into a full LIF network displays 
SSA dynamics with a spectral gap, as predicted by our 
simple rate model, and the paired excitatory/inhibitory 
neurons act as a functional circuit displaying synchronous 
firing behavior. Importantly, note that, in contrast to the 
excitatory clustered scenario, there are no groups with 
dense reciprocal couplings in this network topology. Our 
numerics also show that varying the strength of the func¬ 
tional grouping (Wei = Wje) leads to the emergence of 
SSA dynamics, linked to the presence of the spectral gap 
AA (Fig. 8 A), and to the alignment of the leading Schur 
vectors with the ‘correct’ functional circuits, i.e., pairs of 
excitatory and inhibitory groups together (Figure 8 B). 

It is interesting to note that in this topology there is 
also a group of eigenvalues bounded away in the negative 
direction (see Figure 7D). These modes are the quickest 
decaying and correspond to ‘anti-correlated’ firing states, 
in which excitatory neurons act in synchrony with the in¬ 
hibitory groups to which they are not functionally associ¬ 
ated. Such fast decay reinforces the survival of synchrony 
within the functional groups in the network. Interest¬ 
ingly, these modes were already present in our linear rate 
model: u 4 in (13) shows the same firing pattern with the 
fastest eigenvalue A 4 = —Vk(s — e). 
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Figure 8. The presence of SSA dynamics in LIF networks 
with excitatory to inhibitory feedback loops. A The spike- 
rate variability, which measures the intensity of SSA dynam¬ 
ics, as a function of Wei = Wie (left; dots: raw data from 
simulations, line: mean, shading: standard deviation) and AA 
(right). See text and Fig. 4 for further details. B The first 
principal angle between the subspaces of the principal firing 
patterns and the dominant Schur vectors of the weight ma¬ 
trix W show high alignment. See text and Fig. 5 for further 
details. 


D. SSA dynamics in LIF networks with alternative 
connectivity topologies 

Although perhaps the most intuitive way of generat¬ 
ing SSA dynamics follows from clustering the connec¬ 
tivity of the LIF network, our analysis above suggests 
that alternative types of structural organization support 
SSA dynamics in LIF networks, as long as they are char¬ 
acterized by a gap in their eigenvalue spectrum and a 
measure of block-localization of the Schur vectors. More 
generally, one could conjecture that any low rank pertur¬ 
bation of the weight matrix of a balanced network which 
leads to an eigenvalue gap might be a valid candidate 
for a mechanism to generate SSA in neuronal dynamics. 
A few such network architectures are worth highlighting, 
as they have been considered of particular relevance in a 
neuroscience context. 


1. SSA dynamics in networks with small-world organization 

The first noteworthy example is the broad class of 
small-world like networks [41], since small-world orga¬ 
nizations have been observed in many structural and 
functional neurophysiological networks [42]. While many 
modular networks can have the small-world property [43], 
here we focus on the archetypal small-world structure a 
la Watts and Strogatz [41], which does not display a dis¬ 
tinct modular organization but instead may be seen as 
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Figure 9. Small-world networks can exhibit SSA activity. 
A Small-world network and the resulting raster plot of its 
LIF network simulation. SSA activity is observed, yet with 
less distinct boundaries when compared to the clustered case. 
B Spectrum of the weight matrix of the small-world network 
showing a small eigenvalue gap. C The principal angle be¬ 
tween the dominant subspace spanned by principal compo¬ 
nents of the firing rates and the dominant Schur vectors of 
the weight matrix shows that, as the modularity of the small- 
world rewiring becomes larger, there is alignment between the 
observed dynamics and the slow directions of the weight ma¬ 
trix (dots: raw data from simulations; line: mean; shading: 
standard deviation). 


works possess distinctive spectral properties which have 
important implications for dynamical processes, such as 
synchronization [44]. As Figure 9A shows, small-world 
like LIF networks (see Materials and Methods for details 
of the construction) can indeed display SSA dynamics. 
In line with our arguments above, the weight matrix has 
a leading group of eigenvalues separated from the bulk, 
which dictate the slow dynamics, although in this case 
the gap is not as large and the slow eigenvalues are not 
as tightly bunched (Figure 9B). These separated eigen¬ 
values correspond to a subspace spanned by eigenmodes 
with a slow spatial variation along the ‘backbone’ ring, 
and which lead to the localization and subsequent switch¬ 
ing of the neural activity between subgroups of neurons. 
The principal angle between the firing patterns and the 
dominant Schur vectors (Figure 9C) indicates that the 
firing patterns are again dictated by the spectral proper¬ 
ties of the underlying small-world topology. Note, how¬ 
ever, that due to the lack of hard boundaries in the spec¬ 
tral groupings, both for the eigenvalue structure and the 
smoothly-varying dominant Schur vectors, the SSA dy¬ 
namics (Figure 9A), although present, is not as distinctly 
localized as in the examples above. This behavior is ef¬ 
fectively a result of the heavy overlap induced by the 
underlying lattice-like pristine world in the small-world 
construction [44]. 


2. Lack of SSA dynamics in scale-free networks 

Another class of networks that has attracted tremen¬ 
dous interest is so-called scale-free networks [45]. These 
networks are characterized by a fat-tailed degree distri¬ 
bution, a fact that has been observed empirically for a 
variety of networks. As shown for undirected random 
graphs, the spectrum of scale-free networks may have a 
spectral gap, due to the presence of hubs with very large 
degrees [27]. However, the associated eigenvectors are 
usually localized around the hubs, and as such are not 
able to trigger a consistent pattern of localized activity 
across a group of nodes. In our simulations of scale-free 
LIF networks, we did not observe SSA dynamics: the 
dynamics was essentially concentrated around the hub 
with no transitions in time (z.e., only a core of neurons 
around the hub had consistently elevated firing rates), 
hence no SSA. Indeed, the networks had only one large 
eigenvalue separated from the bulk (due to the high de¬ 
gree of the main hub) and its associated Schur vector was 
highly localized around the hub [27]. An example of this 
hub-centric dynamics can be found in Figure S2 of the 
Supplementary Information. 


3. SSA dynamics with multiple time scales in LIF networks 


a mixture of a fc-nearest neighbor ring lattice and a ran¬ 
dom graph. It is well known that such small-world net- 


It is also possible to enforce hierarchical arrangements 
of the functional wiring mechanisms discussed so far, thus 
allowing for a variety of combinatorial arrangements. De- 
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pending on the number of hierarchical layers, this enables 
the introduction of multiple time scales in the spatio- 
temporal segregated activity. As an illustrative example, 
Figure 10 shows the dynamics of a LIF network with a 
two-level hierarchy of clustered excitatory neurons. This 
network exhibits SSA with two ‘slow’ time scales: a very 
slow switching between the groups of the top level of 
the hierarchy, and the not-so-slow pulsation of activity 
between the subgroups within each of the large groups, 
thus reflecting the second level of the hierarchy. We re¬ 
mark that one could have employed a combination of dif¬ 
ferent ‘functional circuits’ in the individual layers of the 
hierarchy leading to potentially different behaviors. As 
real neural networks display multiple layers of hierarchi¬ 
cal organization [15-18, 46], understanding such schemes 
is of great importance for neural computation, and will 
be addressed in future work. 


III. DISCUSSION 

Answering the question of how the wiring of neural 
circuits governs neural dynamics is a key step in under¬ 
standing neuronal computations. Here, we showed how 
knowledge of certain spectral features of the matrix of 
neuronal connectivity can help understand what dynam¬ 
ics a network can support. 

In this paper we focused on LIF networks exhibiting 
slow-switching assembly (SSA) dynamics, where groups 
of neurons show a sustained and coherent increase in their 
firing rates, with slow switching between epochs of lo¬ 
calized firing in different groups across the network. We 
found that the presence of the slow switching time scale is 
reflected in the spectral properties of the synaptic weight 
matrix: a gap separating the leading eigenvalues together 
with a block-localization of the associated Schur vectors 
on groups of neurons is a key indicator of the presence 
of SSA dynamics in the network. In line with this ob¬ 
servation, multiple gaps in the eigenvalue spectrum are 
indicative of further time scales in the network dynamics 
(see Figure 10). Moreover, when the leading eigenvalue 
becomes larger than one, the SSA dynamics becomes lo¬ 
calized on one of the assemblies (see Figure 6). 

First, we revisited the case of balanced LIF networks 
with clustered excitatory neurons [22] and observed that 
only when there is an eigenvalue gap and Schur block- 
localization does the network display SSA dynamics. 
Further analytical understanding from a stylized firing- 
rate model allowed us to determine that the clustering 
strength drives the development of the eigenvalue gap 
responsible for the slow switching between localized fir¬ 
ing modes. We remark that clustered excitatory connec¬ 
tivity leads to a Hebbian amplification regime [24]: the 
more stable the pattern formation, the slower the dy¬ 
namics of any banded activity. Fast switching between 
well-defined, stable patterns is thus only achievable for 
strong input changes. 

As suggested from our spectral characterization, and 
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Figure 10. SSA activity in hierarchical LIF networks. 
A Schematic of a hierarchically clustered LIF network. Each 
excitatory group consists of two subgroups. B A hierarchical 
modular LIF network with N — 2000 neurons and the result¬ 
ing raster plot of its simulation. SSA activity is observed at 
two time scales, corresponding to the two hierarchical levels 
embedded in the network structure: slow switching between 
the large groups and faster switching between the inner sub¬ 
groups (see inset). C The spectrum of the weight matrix of 
this hierarchical network exhibits two eigenvalue gaps corre¬ 
sponding to the two slow time-scales the network can support. 


confirmed by simulations, we showed that SSA dynamics 
can be achieved through the structured tuning of synap¬ 
tic strengths, rather than through the clustered rewiring 
of the connections. While the equivalence between topo¬ 
logical and weight organizations is clear for a linear sys- 
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tem, such a direct correspondence is not guaranteed a 
priori for a non-linear system. In fact, the clustering of 
weights appears to have a slightly larger influence com¬ 
pared to topological organization in our simulations of 
non-linear LIF systems. More importantly, this dynami¬ 
cal equivalence between clustered topology and clustered 
synaptic weights has potential ramifications for learning 
and synaptic plasticity; it indicates that the alteration of 
the weight structure can lead to the emergence of grouped 
activity without the need for structural rewiring of phys¬ 
ical synaptic connections, thus suggesting a cost-effective 
adaptation to stochastically encode different firing pat¬ 
terns [47]. 

Within LIF networks with clustered excitatory con¬ 
nectivity [22], inhibitory neurons only play a balancing 
background role, whereas experimental studies have re¬ 
vealed a vast diversity of interneuron subtypes that play 
key roles in computation [48, 49]. We thus investigated 
mechanisms in which inhibitory neurons have an active, 
functional role in generating SSA dynamics. We demon¬ 
strated both analytically and via numerical simulations 
how such functional circuits can be constructed with 
the help of a feedback loop between excitatory and in¬ 
hibitory neurons, so that inhibitory neurons themselves 
exhibit SSA dynamics and become an integral functional 
part of the dynamics. As fewer inhibitory neurons are 
present, adapting the weights according to this functional 
co-clustering scheme may also provide a cost effective al¬ 
ternative to generate SSA activity in neuronal networks. 

The emergence of SSA dynamics can be achieved not 
only by clustering excitatory neurons but also by shap¬ 
ing the network structure in several ways. For instance, 
a small-world network topology of the excitatory neurons 
is also able to support SSA dynamics due to the spectral 
properties of small-worlds, which induce a separation of 
the slow eigenvalues and the localization and switching 
of firing activity associated with slowly spatially-varying 
dominant Schur vectors. In contrast, scale-free networks 
lack the block-localization of the Schur vectors on multi¬ 
ple groups which is necessary to consistently induce SSA 
dynamics in LIF networks. Further topologies will be 
explored in the future. In particular, inhibitory neurons 
that inhibit other inhibitory neurons (so-called disinhibi- 
tion patterns [48, 50-54]) provide an interesting example. 
A different avenue might be provided by balanced ampli¬ 
fication mechanisms [24], which could be potentially used 
for creating grouped activity, yet without the introduc¬ 
tion of a slow time scale. Other interesting questions in 
this respect are: which wiring mechanisms provide the 
most economical, or evolutionarily fit, variant [47, 55] to 
induce a given dynamics; and how does the observed di¬ 
versity of interneurons (and the multiple roles they can 
play) relate to this dynamical picture. 

Our work emphasizes the importance of the spectral 
characterization of the weight matrix for the dynam¬ 
ics taking place on these topologies. Although spectral 
properties are also key in characterizing the dynamical 
response of networks operating at criticality [56], our ob¬ 


servations here correspond to a different phenomenon. 
The emergence of a dominant assembly when A max > 1 
leads to a saturation of the network and a decrease of its 
dynamical heterogeneity, in contrast to systems of cou¬ 
pled non-modular excitatory units at criticality, which 
maximize their dynamic range for A max = 1 [56, 57]. 

Our work links up with experimental findings that cor¬ 
tical states can arise spontaneously and exhibit switch¬ 
ing behavior. Recordings from anesthetized cat visual 
cortex have revealed activity states that dynamically 
switch, and these states match closely to recorded ori¬ 
entation maps [58]. Such cortical states likely arise due 
to similarly tuned neurons having higher intra-cortical 
connectivity [11, 12], as is also predicted through mod¬ 
els of orientation maps [59-61]. These observations are 
similar to our theoretical outcomes which indicate that 
dynamic transitioning between different network states 
can be driven and sustained based only on the underly¬ 
ing network topology. Other experiments have identified 
states in CA3 network activity of the hippocampus [62] 
where active cell assemblies exist for tens of seconds be¬ 
fore sharply transitioning to a new state. However, these 
were metastable states (hence unlikely to be reactivated) 
and included a core population of cells consistently active 
in all states. While such complex dynamics are beyond 
the simple models presented here, it may be feasible to 
model such metastable dynamics through more elaborate 
network topologies, which may include synaptic plastic¬ 
ity [63]. Finally, let us remark that our reduced model 
descriptions effectively focused on networks implement¬ 
ing a rate-based coding mechanism and do not include 
spike-timing of cell assemblies, which has been identified 
to be an important component in neuronal computation, 
e.g. in rodents [64, 65], monkeys [66], songbirds [67], and 
grasshoppers [68]. Whether our results can be translated 
to a time-coding regime would be an interesting question 
for future work. 

As connectomics continues to advance the mapping of 
connections and synaptic strengths in wide areas of the 
brain, experimentally obtained weight matrices can be 
analyzed spectrally as described above to determine if 
SSA dynamics can be supported by the networks under 
study. While we focused here on the simplest mecha¬ 
nisms underlying SSA activity, future work could also 
consider the construction of biophysically realistic models 
that take into account the growing literature on the dis¬ 
tribution of synaptic strengths, synaptic contacts, firing 
rates, and other relevant cortical parameters that tend 
to show a lognormal distribution [7, 69-71]. 

Although knowledge of the network structure (connec¬ 
tomics) is not sufficient to predict the dynamics a net¬ 
work circuit will display (for instance, the network dy¬ 
namics can be dominated completely by a strong input), 
it can still give valuable insights on the firing patterns 
the network can support. Conversely, the observation of 
neuronal dynamics alone may not be sufficient to under¬ 
stand neural computation in detail, since different net¬ 
work topologies can yield similar dynamics. Hence, our 
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work hints at how connectomics and neuronal dynamics 
data can provide complementary and intertwined routes 
for systems neuroscientists to study the computational 
principles implemented by the brain. 

Finally, it is important to remark that in order to get 
a fuller picture of the relation between structure and dy¬ 
namical network properties many other factors not con¬ 
sidered here will be of interest, including the distribution 
of inputs, the precise location of synapses on the post- 
synaptic neuron, and the plasticity rules that govern the 
evolution of the network over time [72]. Linking the in¬ 
sight gained from our work to experimental observations 
or to highly detailed computational models [3, 73] would 
thus be a fruitful next step in order to bridge the gap 
between neuronal structure, dynamics, and ultimately 
function. 

IV. MATERIAL AND METHODS 
A. Leaky integrate-and-fire networks 

We simulated leaky-integrate-and-fire (LIF) networks 
where the non-dimensionalized membrane potential of 
each neuron (V*(t), i = 1 ,..., N ) was modeled by: 

= —(Mi - Vi(t)) + (14) 

CLX T rn 

3 

with a firing threshold of 1 and a reset potential of 0. The 
constant input terms fa were chosen uniformly in the in¬ 
terval [1.1,1.2] for excitatory neurons, and in the interval 
[1,1.05] for inhibitory neurons. The membrane time con¬ 
stants for excitatory and inhibitory neurons were set to 
r m = 15 ms and r m = 10 ms, respectively. The refrac¬ 
tory period was fixed at 5 ms for both excitatory and in¬ 
hibitory neurons. Note that although the constant input 
term is supra-threshold, balanced inputs guaranteed that 
the average membrane potential is sub-threshold [22]. 

The network dynamics is captured by the sum in (14), 
which describes the input to neuron i from all other 
neurons in the network. The topology of the network 
is encoded by the weight matrix IF, i.e., Wij denotes 
the weight of the connection from neuron j to neuron 
i, where Wij is zero if there is no connection. Synaptic 

inputs are modeled by g^ 1 (t), which is increased step¬ 
wise instantaneously after a presynaptic spike of neuron 
j (gf /T gf^ 1 + 1) and then decays exponentially ac¬ 
cording to: 

dgf^ 1 e/i 

TeI 1 Wt = ~ 9 i F’ ( 15 ) 

with time constants te = 3 ms for an excitatory inter¬ 
action, and rj = 2 ms if the presynaptic neuron is in¬ 
hibitory. 

Equation (14) can be rewritten in matrix notation as: 

d ^P- = T~ 1 [n-V]+Wg(t), (16) 


where T = diag(r^) and V = [Vi,..., Vn] t and the vec¬ 
tors V, /j,, and g are N x 1 vectors. The spectral analy¬ 
ses in the main text (eigenvalues and Schur vectors) refer 
to the weight matrix IF. Different network topologies 
correspond to different weight matrices W , as explained 
below. 

In all simulations, the ratio of excitatory to inhibitory 
neurons was fixed to be 4 : 1 . Unless otherwise stated, the 
networks comprise N = 2000 units (1600 excitatory, 400 
inhibitory) and were simulated over 20 seconds to calcu¬ 
late the statistics reported. The time step for all simula¬ 
tions was 0.1 ms. All simulations were run in MATLAB 
(version 2011 b or later). 


B. Network Topologies and Weight matrices 

Different network topologies were simulated using dif¬ 
ferent IF matrices but maintaining a general balanced 
network. If not indicated differently, all parameters be¬ 
low correspond to the case N = 2000. 


1. Networks with unclustered balanced connections 

Unless otherwise stated, unclustered connections be¬ 
tween neurons were drawn uniformly at random accord¬ 
ing to the probabilities p EE = 0.2 and p EI = p IE = 
p 11 = 0.5, where the first superscript denotes the desti¬ 
nation and the second superscript denotes the origin of 
the synaptic connection, and E, I stand for an excita¬ 
tory or inhibitory neuron, respectively. Synaptic weights 
between inhibitory neurons had always the weight w 11 = 
—0.0297. The other weight parameters were set to w EI = 
—0.0297, w IE = 0.0074 and w EE = 0.0156, except where 
indicated differently in the scenarios described below. 


2. Networks with clustered excitatory-to-excitatory 
connections 

First, LIF networks with clustered excitatory neurons 
were constructed according to the protocol in Ref. [ 22 ]. 
The connection probabilities between excitatory neurons 
were changed so that neurons had more connections 
within its group than to neurons outside the group. The 
average number of connections was kept the same as for 
an unclustered balanced network ( p EE = 0 . 2 ) while vary¬ 
ing the ratio R E e = Pin/Pout: where p in and p out are the 
probabilities of connectivity within a group and between 
groups, respectively. Ree = 1 corresponds to the un¬ 
clustered balanced network. 

Second, we generated LIF networks where the con¬ 
nectivity between all excitatory neurons was uniform 
( p EE = 0 . 2 , Ree = 1) but the weights were varied by 
changing the ratio Wee = w f n E / w oiE j where wf n E refers 
to synaptic weights within the same group and w E ^ 
refers to synaptic weights to neurons in other groups. 
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The average synaptic weight between excitatory neurons 
was kept at w EE = 0.0156, as for the unclustered case. 
Strictly speaking, this scheme does not correspond to a 
clustering in a topological sense, but rather to a clus¬ 
tering of synaptic weights between excitatory neurons. 
However, in terms of the dynamics, both E-E clustered 
network variants (adjusting probabilities or weights) lead 
to effectively the same results, as discussed in the text. 
Unless otherwise stated, the excitatory neurons were di¬ 
vided into 20 groups of 80 neurons each. 

As the weights in balanced networks should scale 
with 1/y/N [ 22 , 23], when varying the network size (see 
Fig. 2 A) the parameter settings were kept the same but 
all weights were scaled accordingly, e.g., for N = 1000 
the weights were multiplied by y/2 compared to those of 
the network with N = 2000. 


3. Networks with excitatory-to-inhibitory feedback loop 

To construct LIF networks with co-clustered excitatory 
and inhibitory neurons, we kept uniform excitatory-to- 
excitatory and inhibitory-to-inhibitory couplings while 
varying the ratios W IE = /w 1 ^, W E i = 

( w fn/ w out)~ l -> or analogous ratios of connections 
probabilities Ri E ,R E j. The subscript ‘in’ indicates con¬ 
nections within the excitatory/inhibitory pair. Note that 
W E i,R E i are defined by the inverse in/out ratio, since 
they describe an inhibitory effect. To modulate the co¬ 
clustered excitatory-to-inhibitory network dynamics, we 
vary the ratios Wi E ,W E i , while keeping the average 
weights constant. For instance, in the example shown 
in Figure 8 , we kept R E j = Rj E = 2 while simultane¬ 
ously varying the ratios Wj E = W E i between 1 and 5. 
For all networks, we divided the network into 20 groups 
with 80 excitatory and 20 inhibitory neurons each. 


4- Networks with small-world connectivity 

Networks with small-world connectivity between exci¬ 
tatory neurons were constructed as follows. Excitatory 
neurons (N = 1600) were connected within a 40-nearest 
neighbor ring with probability pf^, i.e., neurons i was 
connected to neuron j with probability p EE if |i—j| < 40 
subject to periodic boundary conditions. The connec¬ 
tion probability outside this neighborhood was set to 
p EE t . As in previous models, the average connectivity 
was kept constant at p EE = 0.2 while varying the ra- 
tio .Rfe = Pin /Pout- Hence, by increasing #§£ we 
increase the amount of ‘backbone’ connectivity through 
the (stochastic) 40 nearest-neighbor cycle, while decreas¬ 
ing the number of connections to neurons elsewhere in 
the network. This construction may thus also be inter¬ 
preted as an overlapping clustering and is effectively a 
variant of the small-world scheme introduced by Watts 
and Strogatz [41]. 


5. Networks with scale-free connectivity 

LIF networks with a scale-free topology for the exci¬ 
tatory connections were created by a preferential attach¬ 
ment process [45] using algorithm 5 in [74] with parame¬ 
ter d = 64. The weight matrix was created iteratively by 
adding one neuron at a time, such that the probability 
that it is connected to an existing node is proportional 
to the degree of that node (at the current state), i.e., the 
earlier a node is added, the larger its final degree will be. 
The resulting degree of the network tends to a power-law 
distribution [45]. The value d = 64 was chosen to avoid 
saturation of the dynamics of the network, and this leads 
to an excitatory-excitatory matrix with 35% of the edges 
of our other examples. 

6. Network with hierarchical excitatory connections 

We simulated LIF networks where excitatory neurons 
belonged to a hierarchy of groups (Figure 10) by divid¬ 
ing the excitatory neurons into nested clusters. We im¬ 
plemented two layers in the hierarchy: the excitatory 
neurons were divided into 16 groups, and each of these 
groups was then subdivided into 2 more groups to give a 
total of 32 subgroups. This was achieved by varying the 

ratios R%* = R*£% = P s E u b g roup/pfrou P while 

keeping the average connectivity between excitatory neu¬ 
rons constant. The example in Figure 10 corresponds to 
R f EE = 1.45 and R S EE = 3.7. The weight within the sub¬ 
groups was set at = 0.0163. All other connections 
remained unchanged. 

C. Quantifying SSA dynamics from spike-train LIF 
simulations 

To evaluate the extent to which a network displays 
slow-switching assembly dynamics, we used the following 
two spike-rate variability measures. 

Given a partition of the neurons into c groups, we com¬ 
pute the average spiking frequency fi(t) of the neurons in 
each cluster i E {1,..., c} over non-overlapping windows 
of 100ms. For each time window, we obtain the firing 
rate vector f (t) = ..., / c (t)] T of which we compute 

the standard deviation cr(t). The standard deviations are 
then averaged over the duration of the simulation: 

s= ^E ( 7 W’ ( 17 ) 

t =i 

where T is the total number of time windows in the sim¬ 
ulation. We then obtain a boot-strapped expectation 
(S'shuff) computed by reshuffling neurons at random into 
groups of the same sizes as those in the partition. The 
spike rate variability score is then defined as: 

S = S-(S shufi ), (18) 
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where the average (S' shuff ) is obtained over 10 random 
reshufflings of the neurons. If the network fires uni¬ 
formly (with no localized patterns), S is low; whereas S 
increases if the network displays heterogeneous activity 
aligned with the partition under investigation. 

As explained in the text, we introduced a second spike- 
train variability measure to quantify the variations of the 
group firing-rates across time. This complementary mea¬ 
sure allows us to discern scenarios in which there is a 
group of neurons dominating the firing (thus leading to a 
large variation across groups), but no switching between 
groups. To quantify these effects we use an analogous 
measure to S above. 

Given a partition of the neurons into c groups, we com¬ 
pute the average spiking frequency fi(t) of the neurons in 
each cluster i G {1,..., c} over non-overlapping windows 
of 100ms. For each group i, we compute the standard de¬ 
viation (iT{i) of the vector of coarse-grained firing rates 
across time fi = [/i(ti),..., /i(tjv)]? where tk stands for 
the kth time bin and i = 1,..., c. We then average over 
groups to obtain: 

St = ~ gtW- (19) 

i= 1 

As above, we obtain a boot-strapped expectation (Sj- n huff ) 
by reshuffling neurons at random into groups of the same 
sizes as those in the partition. The spike rate variability 
score over time is then defined as: 

S T = S T - (20) 


for all cases where SSA dynamics is present, both mea¬ 
sures S and St behave consistently. On the other hand, 
St detects the end of the SSA region when, through in¬ 
creased clustering, the dynamics gets localized on one cell 
assembly, as shown in Figure 6 . 


D. Measuring alignment of LIF dynamics with the 

Schur vectors of the weight matrix: the principal 
angle 

First, we find the dominant firing patterns in the net¬ 
work dynamics. We perform simulations of the LIF net¬ 
work and obtain the firing rates of every neuron in 250 ms 
bins to generate aniVxT matrix, where N is the number 
of neurons and T is the number of bins. On this matrix, 
we perform a principal component analysis (PCA) and 
select the first c— 1 principal components {pijqZj 1 2 3 4 5 6 7 8 9 . This 
set of TV-dimensional vectors captures most of the vari¬ 
ability observed in the simulated dynamics. 

We then assess how aligned the c — 1 principal compo¬ 
nents are with the dominant Schur vectors of the weight 
matrix of the network. More precisely, we compute the 
Schur decomposition of the weight matrix W; keep the 
c—1 dominant Schur vectors {iii}?”* associated with the 
eigenvalues with largest real part; and then compute the 
(first) principal or canonical angle 0 between the sub¬ 
spaces spanned by the two sets of vectors V = span{pi} 
and U = spanjui} [38, 39]: 


where the average (5fi luff ) is obtained over 10 random 
reshufflings of the neurons. Again, if the network fires 
homogeneously (with no localized patterns in time), St 
is low; whereas St increases if the network displays het¬ 
erogeneous firing rates over time. 

In Figure SI, we show the behavior of St for the ex¬ 
amples discussed in the main manuscript. As expected, 


cos($) = max 


T 

tr p 


u 


u eU p g V 


( 21 ) 


The first principal angle measures how ‘close’ the ob¬ 
served firing patterns are to the dominant modes com¬ 
puted solely from the weight matrix. If cos (0) « 1, there 
is a large alignment between the span of both subspaces. 
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Appendix A: Supplementary Information - 
Emergence of slow-switching assemblies in 
structured neuronal networks 
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A “ 1000 neurons; 10 groups 
— 2000 neurons; 20 groups 



C Ree 




Figure SI. Spike rate variability across time, as a function 
of topological changes of the network. A Clustering Strength 
Ree vs spike rate variability across time St for varying net¬ 
work sizes with a fixed assembly size. Compare to Fig. 2A 
of the main text. B Clustering Strength Ree vs spike rate 
variability across time St for a fixed network size with vary¬ 
ing assembly sizes. Compare to Fig. 2B of the main text. 
C Weight Clustering Wee vs spike rate variability across time 
St for a LIF network with 2000 neurons (20 groups; topolog¬ 
ical clustering Ree — 1)- Compare to Fig. 4 of the main 
text. D Spike rate variability across time St as a function of 
Wei — Wie for a network with excitatory to inhibitory feed¬ 
back loops. Compare to Fig. 8 of the main text. In all plots 
dots denote raw data from simulations; line: mean, shading: 
standard deviation. 
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Figure S2. Network dynamics for a scale-free network. 
A Raster plot for an example scale-free network. Excita¬ 
tory neurons are ordered according to their (expected) degree: 
high (bottom) to low (top). Note that the spiking activity is 
largely concentrated around the hub (which due to assorta- 
tivity is connected to other nodes of high degree), and there 
are no distinguishable groups or switching behavior. B Spec¬ 
trum of the simulated scale-free network. There is only one 
eigenvalue separated from the main bulk of the spectrum due 
to the presence of the large hub in the network. 












